Parameter mismatch estimation using large deviations from synchronization 



Jupiter Bagaipo 1:2, Q and Juan G. Restrepo 2,3 

1 Department of Physics, University of Maryland, College Park, Maryland 20742, USA 
Department of Mathematics, University of Maryland, College Park, Maryland 20742, USA 
' Institute for Research in Electronics and Applied Physics, 
University of Maryland, College Park, Maryland 20742, USA 
(Dated: February 8, 2008) 

We present a method to determine the relative parameter mismatch in a collection of nearly iden- 
tical chaotic oscillators by measuring large deviations from the synchronized state. We demonstrate 
our method with an ensemble of slightly different circle maps. We discuss how to apply our method 
when there is noise, and show an example where the noise intensity is comparable to the mismatch. 
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The study of networks of coupled dynamical systems 
is an important area of research with applications in di- 
verse fields, ranging from biology to laser physics 0-0 • 
The synchronization of coupled oscillators has been un- 
der extensive study in recent years and, in particular, 
the synchronization of identical oscillators has received 
considerable interest 0, Q . Since it is impossible in prac- 
tice to obtain identical oscillators, the effect of the dif- 
ference in the parameters of the oscillators, or parameter 
mismatch, might be relevant in some applications. It 
might be desired to have dynamical units as similar to 
each other as possible, or to know the characteristics of 
the parameter mismatch in a collection of nearly iden- 
tical systems. In this paper we propose a method to 
use deviations from synchronization to extract informa- 
tion on the parameter mismatch of the coupled dynamical 
units. Existing methods for parameter estimation (see, 
for example, Refs. usually rely on knowledge of 

the typically small synchronization error. Our method 
depends on relatively large deviations from the synchro- 
nized state, and might be useful in cases in which the 
small synchronization error can not be measured accu- 
rately. 

When a number of identical systems are appropriately 
coupled in a network, a solution exists in which the state 
of all oscillators at all times is the same. This is referred 
to as identical synchronization. This concept is useful 
only when the systems are identical. We will deal with 
systems that are nearly, but not exactly, identical. We 
will refer to a situation in which the states of the systems 
are very close to each other as nearly identical synchro- 
nization. A method to determine the stability of the 
synchronous state when the systems are identical, the 
master stability function, has been proposed by Pecora 
and Carroll |9(. 

In the case of nearly identical chaotic systems, the 
nearly synchronized state might be interrupted by rela- 
tively short periods of desynchronization {desynchroniza- 
tion bursts). These bursts develop with spatial patterns 
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on the network. These spatial patterns, and the param- 
eters for which they can be expected, can be predicted 
from the Laplacian matrix describing the network con- 
nections, the master stability function of the attractor, 
and the unstable periodic orbits embedded in it [lOj . The 
spatial patterns of the bursts depend on the parameter 
mismatch of the different systems. We use this fact to 
infer the relative deviation of the parameters of the indi- 
vidual units with respect to their mean from the desyn- 
chronization bursts. The proposed method is as follows. 
The oscillators are connected in such a way that the pa- 
rameter mismatch determines the spatial patterns of the 
desynchronization bursts. As we will see later, one such 
way is all-to-all coupling. The system is set up in a pa- 
rameter region in which desynchronization bursts are ex- 
pected. While a burst is developing, measurements are 
taken of the deviations of the different systems from the 
synchronous state. From these observations, the relative 
deviations of the parameters from the mean are deduced. 



Some limitations of this method are the following. It 
is assumed that measurements can be taken with enough 
precision such that the deviations from the synchronous 
state can be measured in the linear regime. Also, the 
presence of noise affects the spatial patterns of the bursts. 
Although we will describe how to deal with the noise, the 
effectiveness of the method decreases as the ratio of noise 
to mismatch increases. It is also assumed that unavoid- 
able small differences in the way in which the systems 
are connected to each other does not introduce a differ- 
ence between the systems which is of the same order of 
magnitude or larger than the parameter mismatch being 
measured. 



In Sec. I we briefly describe the master stability func- 
tion method and its extension to deal with nearly iden- 
tical systems. In Sec. II we present and illustrate our 
method with an example in the case where the noise is 
negligible. In Sec. Ill we discuss how to deal with the 
noise and show an example. In Sec. IV we present our 
conclusions. 
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I. BACKGROUND 

For simplicity, we will use one dimensional maps for the 
chaotic units. The results generalize to other dynamical 
systems 0, . We consider a model system of N dy- 
namical units, each one of which, when isolated, satisfies 
X^ +1 = F(X^,fXi), where X l n is the value of unit i at 
time n and /Zj is a parameter vector for system i. 

The systems, when coupled, are taken to satisfy (e.g., 



B) 



x; l+1 = F(xi^) - g z\J2 G l3 H{X^ 



(1) 



where Z is a function such that Z(Q) = 0, G is a sym- 
metric Laplacian matrix Gij — 0) describing the net- 
work connections, and H is a function independent of i 
and j. [In our examples, we will take Z{x) = sin(27rx).] 
The constant g determines the strength of the coupling. 

If the systems are identical (i.e., //j = fi for all i), 
there is an exactly synchronized solution of Eqs. QJ, 



X,, 



A .. 



X 



N 



whose time evolu- 



tion is the same as the uncoupled dynamics of a single 
unit, s n+ i = F(s n ), where F(s) — F(s,fj,). The sta- 
bility of the synchronized state can be determined from 
the variational equations obtained by considering an in- 
finitesimal perturbation 5 l from the synchronous state, 



X' 



JY 



•5^4-1 = DFis^Si ~ gZ'(Q) G l3 DH(s n )6l (2) 



Let S = [5 1 , 8 2 , . . . , 5 N ], and define the vector rj — 
[ry 1 , i] 2 , . . . , rj N ] by S = r/L T , where L is the orthogo- 
nal matrix whose columns are the corresponding real 
orthonormal eigenvectors of G; GL — LA, A = 
diag(Xi, A2, . . . , Ajv), where A& is the eigenvalue of G for 
eigenvector k. Then Eqs. (J2J are equivalent to 



r, k = [DF(s n ) - gZ'(0)X k DH( Sn )} v k . 



(3) 



The quantity rj k is the weight of the kth eigenvector of 
G in the perturbation 6. The linear stability of each 
'spatial' mode k is determined by the stability of the 
solution of Eq. Q . By introducing a scalar variable a = 
gZ'(0)\k, the set of equations given by Eq. 10 can be 
encapsulated in the single equation, 



r) n+ i = [DF(s n ) - aDH(s n )] r] n . 



(4) 



The master stability function ^(a) is the largest Lya- 
punov exponent for this equation. This function depends 
only on the coupling function H and the dynamics of an 
individual uncoupled element, but not on the network 
connectivity. The network connectivity determines the 
eigenvalues Xk (independent of details of the dynamics 
of the chaotic units). The stability of the synchronized 



state of the network is determined by ^P* = sup fc $(pAfe), 
where ^P* > indicates instability. 

If the systems are slightly different, one gets instead of 
Eqs. (J2J) the equations 



N 
J = l 



= DF(s n )5 l n - 9 Z'(0) > ^ G l3 DH(s)5l + Q\s n ), 



(5) 

where F(X l ) = F(X\ £)f m/N), and Q l {X l ) = 
F(X l , Hi) — F(X l ) represents the effect of the mismatch 
and is assumed to be small. Terms of order Q5 were ne- 
glected. Defining Q = [Q 1 (s„), Q 2 (s„), . . . , Q N (s n )], we 
obtain an equation analogous to Eq. ©, 

rfc +1 = [DF(s n ) - gZ'(0)X k DH(s n )} rfc + (QL)\ (6) 

where (QL) k is the kth element of the vector QL. The 
Lyapunov exponent for the solution of Eq. is hi, = 
'i'(gXk)- Assuming a solution of Eq. |J§} to have the same 
average damping as that for Eq. J3J), then, if h k is nega- 
tive for all modes, the amplitude of rjk can be estimated 
as 



(\(QL) k \) 

1 _ e -\h k \ ' 



(7) 



(For example, if we model Eq. © by the simple sys- 
e~ h f]n + Q, then r] n satisfies, as n — > 00, 
See The largest Lyapunov expo- 



tem r] n+1 

q 

nent hj~ above corresponds to a typical trajectory in the 
chaotic attractor. However, the Lyapunov exponent for 
unstable periodic orbits embedded in the attractor might 
be larger. Assume that one of these periodic orbits has 
a positive Lyapunov exponent and the attractor has a 
negative Lyapunov exponent. In this case, most of the 
time the amplitude of rjk will be very small and given 
approximately by Eq. (JJJ. Eventually the trajectory s n 
will get very close to this transversally unstable periodic 
orbit. While it is close to this orbit, 77 in Eq. © is no 
longer damped and gets exponentially amplified with the 
Lyapunov exponent of the unstable periodic orbit. The 
deviation from the synchronized state becomes large, pro- 
ducing a desynchronization burst. If there are no other 
attractors, the system returns to the synchronized state 
and the process repeats. Desynchronization bursts can 
be expected when the master stability function for a typ- 
ical trajectory is negative for all modes, and there is 
an embedded unstable periodic orbit which has a pos- 
itive master stability function for at least one mode (i.e., 
^(t/Afc) > for some k, where ^ is the master stability 
function of one of the embedded periodic orbits). 

We will now use the fact that modes with the same 
eigenvalue have the same stability. For simplicity, as- 
sume that the coupling is all to all, so that all the modes 
(except the mode in the synchronization manifold, which 
has zero eigenvalue) have the same eigenvalue, Xk = N. 
Eq. J7J implies that the coefficients in the eigenvector de- 
composition of the deviations from synchrony (r] k ) are, on 
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average, proportional to those for the deviations of the 
mismatch parameters from their mean \(QL) k ]. It fol- 
lows that the mismatch vector Q is proportional to the 
vector 5 n while these approximations are valid. If the 
vector S n is measured, the deviations of the mismatch 
from its mean can be determined approximately up to 
an unknown scaling factor. 

For this method to work, the measurements need to 
be made when the system is still in the linear regime. 
Since it is assumed that there is a limitation in the mea- 
surement accuracy, S n needs to be small enough to guar- 
antee linear behavior, but large enough to be measured. 
One can thus set up the system so that desynchroniza- 
tion bursts are expected and make measurements while 
a desynchronization burst is developing. If the system 
allows continuous tuning of the coupling strength, one 
could also increase it so that the synchronous state be- 
comes unstable and take measurements as the system 
desynchronizes. 



II. PARAMETER MISMATCH WITHOUT 
NOISE 

To illustrate our method, we use the circle map, de- 
scribed by the equation 

9 n +i = [0 n + w + Ksin27r6'„] mod 1. (8) 

We choose the parameters to be u = ^ and n — 
These parameters produce a chaotic attractor in 
9 e [0.21,0.47]. We found the embedded periodic orbits 
up to period four. To determine the orbits of period p 
we used Newton's method to find the roots of 9 = f p (9), 
where f(9) is described by Eq. (JSJ) and f p denotes the p 
times composition of /. Eliminating all the orbits outside 
of the attractor, we found one period 1 orbit, two period 
2 orbits, and one period 4 orbit. We show in Fig. 2] 
the master stability functions of the orbits found, and 
the master stability function of the attractor. Here a = 
27rgAfc, where A& is the fcth eigenvalue of the coupling 
matrix and g is the global strength of the coupling. 

For definiteness, we assume a network that is coupled 
all to all. This means that for a network of N systems, 
an element in the coupling matrix Gij is given by 

This matrix has two distinct eigenvalues, Ao = and 
Afc = N for k = 1, 2, . . . , N — 1. We ignore the eigen- 
value since this corresponds to a perturbation in which 
all of the systems are displaced by the same amount (thus 
they remain synchronized). Due to the lack of other dis- 
tinct eigenvalues, it is easy to pick an a that produces 
desynchrozination bursts. For our map the region where 
bubbling is expected is where 0.9 < 2-xgXk < 2.1 (note 
that for our example A& is the same for all k). 



4J(a) 




FIG. 1: Master stability function, ^>{pt), for a typical trajec- 
tory in the attractor (continuous curve), for the period 1 orbit 
(dotted curve), for the period 2 orbit (dashed-dotted curve), 
and for the period 4 orbit (dashed curve). 



We now present an example of the method where noise 
is negligible. The mismatch is chosen to be in k since it 
has a more complicated effect than mismatch in oj. The 
coupled systems can then be described by the general 
equation 

#n+i = & + w + {k + 6k 1 ) sin 2n9l - mod 1, (10) 

where Sn 1 is the mismatch in system i, = 
gsin ^27r J2f=i Gij@hj i an d i,j — 1,2,. . .,N are indices 
representing the ith and jth system in the network [cf. 
Eq. ||TJ]. The term $ determines the coupling of the net- 
work. We chose N = 5 systems, g = 0.0477 as the global 
coupling, and 5k — [4,-1,2,-6,-2] x 10~ 6 as the val- 
ues for our examples. It should, however, be noted that 
this method works for any number of systems (with g 
being adjusted accordingly) and mismatch of any size, 
although, if the mismatch becomes small, the waiting 
time for a desynchronization burst becomes large. The 
waiting time can be adjusted by changing the values of 
N, <?, and Sn. 

Defining A9 z n = sin2n(6 l n - 9 n ), where 9 n = 

jj YliLi n , we plot A9 J n versus n and look for desynchro- 
nization bursts in the network. In Fig. [21 we show the 
time evolution of A9 l n near a desynchronization burst. 
Our interest is in the vector 9i where I is the first 
time that maxi{9 l n } is in the sampling region, defined 
as 0.4 < |A0*J < 0.6 (see Fig. EJ. This region is de- 
termined by the limitations on the ability to accurately 
measure A9 l n and the dynamics of the system considered. 
The latter exists because the master stability function 
method relies on the systems being close to synchroniza- 
tion. During the desynchronization burst, the difference 
in the systems can be so large that the linearization used 
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FIG. 3: Superimposed plot of a(Sn 1 — 5k) (dashed line with 
FIG. 2: Plot of A6 l n vs. n near a region of desynchronization square markers) and 9] — 6i (solid line with circle markers) 
burst. The arrow points to the maximum of A(9^ that is versus i with a ~ 1.5 x 10 4 . 
within the sampling region. 



in the analysis of Sec. I no longer applies. We determined 
that the upper bound to this region in the circle map is 
|A0jJ w 0.6 or \9 l n -9 n \ ~ 0.10. The lower bound was ar- 
bitrarily chosen as representing the accuracy of the mea- 
surements, which we assume is not enough to measure 
the mismatch directly. Generally the method becomes 
more effective the smaller the lower bound is. 

According to the previous section, at time I we should 
have approximately 6) — 9i oc 5k 1 — 5k, where 5k = 
jr YliLi We can then obtain the relative deviations 
of the mismatch parameters, 5k % , by measuring the much 
larger values of 9\ — 9\. In Fig.|3|we show a superimposed 
plot of 8\ — 6i and cl(5k 1 — 5k) versus i where a, the scaling 
factor, minimizes Yli=i[{@i ~ &l) ~ o-{5k 1 — 5k)] 2 . InFig.|3| 
we calculated a » 1.5 x 10 4 and this corresponds to the 
amplihcation of the mismatch. It should be noted that 
the sign of a is undetermined unless we have knowledge 
of 5k 1 . We see from the figure that a(5K l — 5k) w 9\ — 9i. 

The definition of the sampling region is somewhat ar- 
bitrary and it may occur that nonlinear effects still play 
a role in the resulting spatial pattern of the burst. In 
fact, in Fig. |21 we observe that there are still small devi- 
ations from the real mismatch pattern. In order to take 
this into account, we can take the average over various 
bursts. In the next section, we will discuss how to ap- 
propriately take the average, and we will also deal with 
the effects of noise. 



III. PARAMETER MISMATCH WITH NOISE 

After learning from the simpler model in Sec. II, we 
can now analyze a more realistic situation. The method 
proposed and explained in the previous section applies to 



a similar network with noise, but there are a few adjust- 
ments to be made. We use the same model described by 
Eq. Ijl0|l . except we modify it to 

9l +1 = [6£+w+(re+tfK i )sin2 7 r6£-<+4] mod 1, (11) 

where 5 k 1 , Q> l n arc defined in the same way as before, and 
e l n is a random variable uncorrelated at different i and n 
and simulates the noise. In our example, we choose e l n 
uniformly from the interval [— 10 -5 , 10 -5 ] (note that the 
noise and mismatch 5k are of comparable size). 

As mentioned in the previous section, nonlinear ef- 
fects might produce deviations from the simple relation 
6\ — 9i = ci(5k 1 — 5k). We assume that the effects of the 
nonlinearity and the noise can be represented by a ran- 
dom variable a 1 , such that 6\ —6i = a(5K l — 5k) + a 1 . We 
furthermore assume that <r l has zero mean. Under these 
assumptions, the mismatch is given by 




where the brackets represent an average over realizations 
of a . Because of the definition of the sampling region, 
the scaling factors for different samples will have similar 
magnitude, but possibly different sign. We thus get ap- 
proximately, assuming the sign of a is independent of a 1 , 

5k 1 -5k tx (sign(a)(^ - 0;)). (13) 

Since the sign of a is unknown, we use a least-squares op- 
timization to find the signs which minimize the dispersion 
from the mean. More precisely, if we have M samples of 
9{s, we can define an average to be = -h J2m=i ftm9m, 
where {fi m } is a sequence of l's and — l's and 9 m is the 
vector [9\ — 9i] for the mth sample. Note that is an 
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iV-dimensional vector and that its ith component is an 
average of the ith component of the M samples of 0;'s. 
If we then minimize the error, defined by 

1 M 

errors — ^ \\(3 m 6 m - Q\\ 2 , (14) 

m— 1 

we can find an optimal sequence of l's and — l's, which 
we shall call (3*, that ensures most of the Ois are oriented 
the same way. 

To minimize this error and optimize (3 it is impractical 
to test all possibilities since there are 2 M different /3's. 
We can, however, follow an algorithm starting with a 
randomly generated /3. At each iterate, we generate three 
new /3's. The first, f3\, is a new random sequence, /?2 is 
(3 altered such that the signs of 1% of the sequence are 
changed, and (3% is defined in a similar way but with 5% 
of the signs changed. We can then compare /3i, /3 2 , (3 3 , 
and (3 and determine which one has the smaller error 
determined by Eq. I|14fl . The one with the smallest error 
is then redefined as (3 and the process is repeated until 
an approximation to (3* , which we denote as (3, is found. 
We can then define 0. = jgEm = iWm- 
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FIG. 4: Superimposed plot of AiSn 1 — Sk) (dashed line with 
square markers) and Ql — 0» (solid line with circle markers) 
versus i with A ~ — 5 x 10 3 . 



We show in Fig.0]a superimposed plot of 6* — 0* and 
A(8k 1 — 5k) versus i where A minimizes X)i=i[(®* — — 
A(8k 1 — 8n)] 2 . To obtain 9* we repeated the process of 
optimization for M — 1000 10 6 times. According to the 
discussion above, we should have approximately @\ — 
0* oc 8k 1 — Sk where 8* = ^ J2i=i ®*- Indeed we see 
that A(Sk 1 — Sk) ps — 9* even when the noise was 
comparable to the mismatch [e ~ Sk in Eq. i|ll|) ]. 
IV. DISCUSSION 

We have presented a method to use large deviations 
from synchronization in order to determine the character- 
istics of the parameter mismatch in a collection of nearly 
identical chaotic dynamical systems. It has been noted 
that knowledge and manipulation of the mismatch pat- 
terns can be advantageous in order to improve the quality 
of the synchronization |10| . The main advantage of our 
method is that it only requires direct knowledge of the 
synchronization error when it is large enough to be mea- 
sured. Furthermore, in principle, there are no limitations 
on the number of systems it can handle. On the other 
hand, the method only provides the relative deviations 
from the mean and has yet to be extended to systems 
with comparable mismatch in different parameters. 

We have demonstrated our method by determining the 
relative parameter mismatch in an ensemble of 5 circle 
maps. By measuring the large deviations from the syn- 
chronized state that occur during a desynchronization 
burst, we were able to determine the very small relative 
differences in parameters (see Fig.[3J). We considered the 
presence of noise, and dealt with it by suitably averag- 
ing the measurements taken for various desynchroniza- 
tion bursts. For a noise comparable to the mismatch 
we were able to determine the relative parameter mis- 
match by averaging 1000 realizations (see Fig. For 
both situations, we were able to determine the relative 
parameter mismatch from measurable values even when 
the mismatch itself was assumed to be immeasurable. 
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